Defining functions and packages

# 1. Source functions ----

source("functions/graphical_par.R")
source("functions/theme_graph.R")
source("functions/pred_vbgf.R")

# 2. Required packages ----

library(tidyverse)
library(kableExtra)
library(formattable)
library(readxl)
library(plotly)
library(rfishbase)

# 3. Set theme_graph() as the default ggplot theme ----

theme_set(theme_graph())

# 4. Load data ----

data_complete <- read.csv("./../data/02_back-calculated-size-at-age_morat-et-al.csv")

1. Data exploration

1.1 Tables

1.1.1 Data with back-calc.

data_complete %>% 
  filter(!is.na(Li_sp_m)) %>% 
  summarize_at(vars("ID", "Species", "Family"), n_distinct, na.rm = TRUE) %>% 
  bind_rows(data_complete %>%
              filter(!is.na(Li_sploc_m)) %>% 
              summarize_at(vars("ID", "Species", "Family"), n_distinct, na.rm = TRUE), .) %>% 
  bind_rows(data_complete %>%
              summarize_at(vars("ID", "Species", "Family"), n_distinct, na.rm = TRUE), .) %>% 
  mutate(Type = c("All data", "Back-calculated sp. loc.", "Back-calculated sp."), .before = 1) %>% 
  kable(., 
        col.names = c("", "Individual", "Species", "Family"), 
        caption = "Table 1. Comparison of numbers of individuals, species and family for the overall dataset and for back-calculated data (by species, and by species and location).") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"))
Table 1. Comparison of numbers of individuals, species and family for the overall dataset and for back-calculated data (by species, and by species and location).
Individual Species Family
All data 855 51 16
Back-calculated sp. loc. 669 44 16
Back-calculated sp. 710 45 16

1.1.2 Individual by species

data_complete %>%
  filter(!is.na(Li_sp_m)) %>% 
  group_by(Species) %>% 
  summarise(n = length(unique(ID)),
        min_lencap = round(min(Lcpt, na.rm = TRUE), 0),
        max_lencap = round(max(Lcpt, na.rm = TRUE), 0),
        max_age = max(Agei, na.rm = TRUE)) %>% 
  kable(., 
        col.names = c("Species", "n", "Min length (TL, mm)", "Max length (TL, mm)", "Age max (years)"),
        caption = "Table 2. Number of individuals, minimum and maximum total length (mm) and maximum age, by species. Only individuals with back-calculated data by species are shown.") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1, italic = T)
Table 2. Number of individuals, minimum and maximum total length (mm) and maximum age, by species. Only individuals with back-calculated data by species are shown.
Species n Min length (TL, mm) Max length (TL, mm) Age max (years)
Abudefduf sexfasciatus 13 61 172 13
Acanthurus achilles 16 138 246 27
Acanthurus lineatus 8 123 362 23
Acanthurus nigricans 8 150 210 9
Acanthurus pyroferus 6 195 230 19
Acanthurus triostegus 18 147 194 12
Balistapus undulatus 17 70 283 18
Caranx melampygus 7 368 712 15
Centropyge bispinosa 7 36 80 11
Centropyge flavissima 26 48 150 27
Cephalopholis argus 41 120 450 21
Cephalopholis urodeta 8 150 214 17
Chaetodon citrinellus 7 64 101 6
Chaetodon ornatissimus 10 122 175 10
Chlorurus spilurus 31 183 344 16
Chromis iomelas 15 41 61 5
Chromis viridis 6 115 148 9
Ctenochaetus marginatus 15 99 270 13
Ctenochaetus striatus 26 153 218 17
Dascyllus aruanus 20 51 71 7
Dascyllus flavicaudus 8 89 104 13
Epibulus insidiator 17 157 350 16
Epinephelus fasciatus 10 123 256 13
Epinephelus hexagonatus 14 160 268 14
Epinephelus merra 46 108 261 17
Epinephelus polyphekadion 14 252 514 20
Gnathodentex aureolineatus 8 197 255 17
Gymnosarda unicolor 4 499 985 7
Halichoeres trimaculatus 7 123 186 5
Lutjanus fulvus 12 195 289 21
Lutjanus kasmira 10 115 326 30
Mulloidichthys flavolineatus 14 257 311 6
Myripristis berndti 30 147 279 24
Naso lituratus 17 228 432 10
Odonus niger 17 176 350 16
Ostorhinchus angustatus 15 46 83 6
Ostorhinchus apogonoides 27 58 96 7
Parupeneus barberinus 6 228 496 6
Plectropomus laevis 30 399 957 22
Pristiapogon taeniopterus 28 53 104 8
Sargocentron microstoma 22 62 204 13
Scarus psittacus 7 195 311 6
Siganus argenteus 11 236 388 13
Stegastes nigricans 13 55 138 13
Zebrasoma scopas 18 110 166 24

1.1.3 Individuals by location

data_complete %>%
  filter(!is.na(Li_sp_m)) %>% 
  group_by(Species, Location) %>% 
  summarise(n = length(unique(ID))) %>% 
  pivot_wider(names_from = Location, values_from = n) %>% 
  mutate(Gambiers = cell_spec(Gambiers, "html", color = ifelse(is.na(Gambiers), "white", "#446CB3")),
         Hao = cell_spec(Hao, "html", color = ifelse(is.na(Hao), "white", "#446CB3")),
         Marquesas = cell_spec(Marquesas, "html", color = ifelse(is.na(Marquesas), "white", "#446CB3")),
         Moorea = cell_spec(Moorea, "html", color = ifelse(is.na(Moorea), "white", "#446CB3")),
         Manuae = cell_spec(Manuae, "html", color = ifelse(is.na(Manuae), "white", "#446CB3")),
         Tuamotu = cell_spec(Tuamotu, "html", color = ifelse(is.na(Tuamotu), "white", "#446CB3"))) %>% 
  kable(., 
        format = "html", escape = FALSE,
        caption = "Table 3. Number of individuals of each species by location. Only individuals with back-calculated data (by species) are shown.") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1, italic = T)
Table 3. Number of individuals of each species by location. Only individuals with back-calculated data (by species) are shown.
Species Gambiers Moorea Manuae Marquesas Tuamotu Hao
Abudefduf sexfasciatus 11 2 NA NA NA NA
Acanthurus achilles 6 NA 10 NA NA NA
Acanthurus lineatus NA NA NA 8 NA NA
Acanthurus nigricans NA NA NA 8 NA NA
Acanthurus pyroferus NA NA NA 6 NA NA
Acanthurus triostegus 3 11 NA 4 NA NA
Balistapus undulatus NA 12 NA 5 NA NA
Caranx melampygus NA 5 NA NA 2 NA
Centropyge bispinosa NA 7 NA NA NA NA
Centropyge flavissima 9 5 NA 12 NA NA
Cephalopholis argus 6 10 NA NA 10 15
Cephalopholis urodeta 8 NA NA NA NA NA
Chaetodon citrinellus NA 7 NA NA NA NA
Chaetodon ornatissimus 4 NA NA 6 NA NA
Chlorurus spilurus 13 16 NA NA 2 NA
Chromis iomelas NA 15 NA NA NA NA
Chromis viridis 6 NA NA NA NA NA
Ctenochaetus marginatus NA NA NA 15 NA NA
Ctenochaetus striatus 11 15 NA NA NA NA
Dascyllus aruanus NA 20 NA NA NA NA
Dascyllus flavicaudus 8 NA NA NA NA NA
Epibulus insidiator 12 5 NA NA NA NA
Epinephelus fasciatus NA NA NA 10 NA NA
Epinephelus hexagonatus NA NA NA NA NA 14
Epinephelus merra 13 20 NA NA NA 13
Epinephelus polyphekadion NA NA NA NA NA 14
Gnathodentex aureolineatus 8 NA NA NA NA NA
Gymnosarda unicolor 4 NA NA NA NA NA
Halichoeres trimaculatus 6 1 NA NA NA NA
Lutjanus fulvus NA 6 NA NA 6 NA
Lutjanus kasmira 2 NA NA 8 NA NA
Mulloidichthys flavolineatus NA 14 NA NA NA NA
Myripristis berndti 7 15 NA 8 NA NA
Naso lituratus 11 6 NA NA NA NA
Odonus niger NA 13 NA 4 NA NA
Ostorhinchus angustatus NA 15 NA NA NA NA
Ostorhinchus apogonoides NA NA NA 27 NA NA
Parupeneus barberinus 5 NA NA 1 NA NA
Plectropomus laevis 12 NA NA NA NA 18
Pristiapogon taeniopterus NA NA NA 28 NA NA
Sargocentron microstoma NA 17 NA 5 NA NA
Scarus psittacus NA 2 NA NA 5 NA
Siganus argenteus 6 2 NA NA 3 NA
Stegastes nigricans 3 10 NA NA NA NA
Zebrasoma scopas 5 13 NA NA NA NA

1.1.4 Number of NA

# 1. Number of NA by variable

data_complete %>%
  summarise_all(~(sum(is.na(.)))) %>% 
  t(.) %>% 
  as.data.frame() %>% 
  kable(., col.names = c("NA"), caption = "Table 4. Number of rows with Non-Available (NA) data by variable") %>% 
  kable_styling(bootstrap_options = c("striped", "hover"))
Table 4. Number of rows with Non-Available (NA) data by variable
NA
Family 0
Genus 0
Species 0
ID 0
Agei 0
Ri 387
Agecpt 0
Rcpt 0
Lcpt 0
L0p 0
R0p 2811
Li_sp_m 410
Li_sp_sd 410
Li_sploc_m 757
Li_sploc_sd 757
Weight 603
Location 0
Observer 0

1.2 Plot

1.2.1 Agei vs Li (sp.)

data_complete %>% 
  filter(!is.na(Li_sp_m)) %>% # Remove individuals with non back-calculated data
  ggplot(., aes(x = Agei, y = Li_sp_m)) +
  geom_point(color = col_color_graph, fill = col_fill_graph, size = 1, shape = 21)+
  theme(strip.text.x = element_text(size = 8, face = "italic"),
        strip.background = element_rect(colour = "black", fill = col_facet, size = 1),
        plot.title = element_text(colour = col_color_graph),
        plot.subtitle = element_text(colour = "black"),
        strip.text.y = element_text(angle = 360, face = "italic")) +
  facet_wrap(~Species, scales = "free", ncol = 5) +
  labs(x = "Age (years)", y = "Length (mm)")
Figure 1. Relation between length (TL, mm) and age, by species.

Figure 1. Relation between length (TL, mm) and age, by species.

1.2.2 Agei vs Li (sp. and loc.)

data_complete %>% 
  filter(!is.na(Li_sploc_m)) %>% # Remove individuals with non back-calculated data
  ggplot(., aes(x = Agei, y = Li_sploc_m)) +
  geom_point(color = col_color_graph, fill = col_fill_graph, size = 1, shape = 21)+
  theme(strip.text.x = element_text(size = 8),
        strip.background = element_rect(colour = "black", fill = col_facet, size = 1),
        plot.title = element_text(colour = col_color_graph),
        plot.subtitle = element_text(colour = "black"),
        strip.text.y = element_text(angle = 360, face = "italic")) +
  facet_grid(Species~Location, scales = "free") +
  labs(x = "Age (years)", y = "Length (mm)")
Figure 2. Relation between length (TL, mm) and age, by species and location.

Figure 2. Relation between length (TL, mm) and age, by species and location.

2. VBGF results

2.1 Growth parameters (sp.)

# 1. Growth parameters (on back-calculated data by species) ----

read.csv("./../data/03_back-calculated_vbgf_predictions_sp.csv") %>% 
  mutate(Estimate = paste0(round(Estimate, 3), " (", round(Est.Error, 3), ")")) %>% 
  pivot_wider("Species", names_from = Parameter, values_from = Estimate) %>% 
  select(Species, linf, k, t0) %>% 
  arrange(Species) %>% 
  kable(., 
        col.names = c("Species", "Linf", "K", "t0"), 
        caption = "Table 5. Values of Von Bertalanffy parameters estimated on back-calculated data through Bayesian framework by species. Standard deviation are the values in parentheses. Linf is expressed in TL and cm.") %>% 
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1, italic = T)
Table 5. Values of Von Bertalanffy parameters estimated on back-calculated data through Bayesian framework by species. Standard deviation are the values in parentheses. Linf is expressed in TL and cm.
Species Linf K t0
Abudefduf sexfasciatus 14.851 (0.394) 1.101 (0.106) -0.026 (0.024)
Acanthurus achilles 18.499 (0.81) 1.185 (0.149) -0.023 (0.024)
Acanthurus lineatus 24.677 (2.411) 1.044 (0.244) -0.029 (0.037)
Acanthurus nigricans 17.114 (0.988) 1.485 (0.229) -0.011 (0.016)
Acanthurus pyroferus 19.594 (1.028) 0.783 (0.137) -0.072 (0.069)
Acanthurus triostegus 15.717 (0.339) 0.796 (0.06) -0.051 (0.028)
Balistapus undulatus 18.987 (1.902) 0.338 (0.059) -0.282 (0.095)
Caranx melampygus 54.386 (4.05) 0.246 (0.039) -0.192 (0.083)
Centropyge bispinosa 5.24 (0.897) 1.099 (4.24) -0.067 (0.036)
Centropyge flavissima 10.036 (0.654) 0.4 (0.067) -0.222 (0.081)
Cephalopholis argus 29.978 (1.244) 0.397 (0.021) -0.251 (0.042)
Cephalopholis urodeta 19.922 (1.553) 0.153 (0.024) -0.569 (0.113)
Chaetodon ornatissimus 14.554 (0.681) 1 (0.141) -0.025 (0.036)
Chlorurus spilurus 23.283 (0.983) 1.172 (0.126) -0.017 (0.019)
Chromis iomelas 5.002 (0.364) 1.292 (0.247) -0.053 (0.021)
Chromis viridis 12.002 (1.096) 0.942 (0.241) -0.05 (0.052)
Ctenochaetus marginatus 20.925 (1.501) 0.424 (0.042) -0.134 (0.044)
Ctenochaetus striatus 17.145 (0.402) 0.849 (0.053) -0.045 (0.022)
Dascyllus aruanus 5.596 (0.2) 0.963 (0.1) -0.05 (0.025)
Dascyllus flavicaudus 8.977 (0.42) 0.499 (0.055) -0.153 (0.059)
Epibulus insidiator 22.532 (1.218) 0.565 (0.061) -0.113 (0.041)
Epinephelus fasciatus 20.424 (1.613) 0.715 (0.109) -0.038 (0.03)
Epinephelus hexagonatus 21.301 (0.928) 0.637 (0.068) -0.067 (0.042)
Epinephelus merra 17.755 (0.759) 0.854 (0.057) -0.045 (0.019)
Epinephelus polyphekadion 38.922 (2.445) 0.281 (0.034) -0.236 (0.074)
Gnathodentex aureolineatus 20.184 (0.928) 0.549 (0.079) -0.154 (0.077)
Halichoeres trimaculatus 17.999 (1.763) 0.645 (0.145) -0.036 (0.058)
Lutjanus fulvus 23.88 (0.928) 0.674 (0.082) -0.086 (0.049)
Mulloidichthys flavolineatus 25.756 (0.53) 1.601 (0.109) -0.008 (0.009)
Myripristis berndti 20.121 (0.761) 0.419 (0.033) -0.235 (0.052)
Naso lituratus 31.222 (1.786) 1.815 (0.237) -0.004 (0.009)
Odonus niger 25.063 (1.219) 0.46 (0.053) -0.141 (0.05)
Ostorhinchus angustatus 7.229 (0.348) 0.681 (0.085) -0.111 (0.036)
Ostorhinchus apogonoides 8.243 (0.223) 0.836 (0.053) -0.072 (0.014)
Plectropomus laevis 71.926 (3.013) 0.317 (0.029) -0.194 (0.056)
Pristiapogon taeniopterus 8.013 (0.272) 0.78 (0.068) -0.093 (0.027)
Sargocentron microstoma 15.714 (0.815) 0.701 (0.077) -0.049 (0.025)
Siganus argenteus 29.483 (2.597) 0.728 (0.143) -0.036 (0.044)
Stegastes nigricans 11.721 (0.543) 0.452 (0.057) -0.16 (0.077)
Zebrasoma scopas 13.319 (0.283) 1.154 (0.089) -0.022 (0.018)

2.2 Growth parameters (sp. loc.)

# 1. Growth parameters (on back-calculated data by species and locations) ----

read.csv("./../data/03_back-calculated_vbgf_predictions_sploc.csv") %>% 
  mutate(Estimate = paste0(round(Estimate, 3), " (", round(Est.Error, 3), ")")) %>% 
  pivot_wider(c("Species", Location), names_from = Parameter, values_from = Estimate) %>% 
  select(Species, Location, linf, k, t0) %>% 
  arrange(Species) %>% 
  kable(., 
        col.names = c("Species", "Location", "Linf", "K", "t0"), 
        caption = "Table 6. Values of Von Bertalanffy parameters estimated on back-calculated data through Bayesian framework by species and location. Standard deviation are the values in parentheses. Linf is expressed in TL and cm.") %>% 
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1, italic = T)
Table 6. Values of Von Bertalanffy parameters estimated on back-calculated data through Bayesian framework by species and location. Standard deviation are the values in parentheses. Linf is expressed in TL and cm.
Species Location Linf K t0
Abudefduf sexfasciatus Gambiers 14.837 (0.393) 1.103 (0.104) -0.027 (0.024)
Acanthurus achilles Manuae 19.44 (0.642) 1.203 (0.14) -0.021 (0.025)
Acanthurus lineatus Marquesas 24.726 (2.413) 1.044 (0.255) -0.029 (0.039)
Acanthurus nigricans Marquesas 17.06 (1.02) 1.495 (0.226) -0.011 (0.015)
Acanthurus pyroferus Marquesas 19.561 (1.02) 0.788 (0.142) -0.07 (0.068)
Acanthurus triostegus Moorea 15.673 (0.476) 0.742 (0.073) -0.067 (0.038)
Balistapus undulatus Marquesas 19.329 (2.174) 0.413 (0.151) -0.254 (0.173)
Balistapus undulatus Moorea 17.235 (1.905) 0.376 (0.066) -0.237 (0.096)
Caranx melampygus Moorea 54.835 (4.393) 0.242 (0.042) -0.167 (0.102)
Centropyge bispinosa Moorea 5.253 (0.849) 0.926 (0.565) -0.066 (0.036)
Centropyge flavissima Marquesas 6.773 (0.573) 1.076 (0.209) -0.033 (0.021)
Centropyge flavissima Moorea 8.175 (0.757) 1.262 (0.329) -0.028 (0.026)
Centropyge flavissima Gambiers 12.21 (1.244) 0.204 (0.041) -0.834 (0.252)
Cephalopholis argus Moorea 22.852 (2.469) 0.728 (0.149) -0.068 (0.052)
Cephalopholis argus Tuamotu 27.447 (1.451) 0.439 (0.064) -0.225 (0.099)
Cephalopholis argus Gambiers 40.416 (1.071) 0.297 (0.031) -0.45 (0.14)
Cephalopholis argus Hao 32.317 (1.775) 0.401 (0.035) -0.184 (0.058)
Cephalopholis urodeta Gambiers 20.981 (1.278) 0.141 (0.021) -0.604 (0.102)
Chaetodon ornatissimus Marquesas 14.599 (1.148) 1.077 (0.261) -0.022 (0.043)
Chlorurus spilurus Moorea 22.998 (1.313) 1.245 (0.214) -0.016 (0.03)
Chlorurus spilurus Gambiers 22.585 (1.059) 0.925 (0.12) -0.037 (0.032)
Chromis iomelas Moorea 4.994 (0.355) 1.301 (0.247) -0.052 (0.021)
Chromis viridis Gambiers 11.953 (1.072) 0.951 (0.228) -0.048 (0.051)
Ctenochaetus marginatus Marquesas 20.934 (1.533) 0.424 (0.043) -0.135 (0.046)
Ctenochaetus striatus Moorea 17.326 (0.566) 1.142 (0.098) -0.021 (0.016)
Ctenochaetus striatus Gambiers 16.99 (0.543) 0.566 (0.051) -0.112 (0.049)
Dascyllus aruanus Moorea 5.589 (0.201) 0.964 (0.101) -0.05 (0.025)
Dascyllus flavicaudus Gambiers 8.992 (0.403) 0.498 (0.054) -0.152 (0.059)
Epibulus insidiator Gambiers 23.859 (1.327) 0.539 (0.066) -0.123 (0.052)
Epinephelus fasciatus Marquesas 20.354 (1.664) 0.721 (0.117) -0.038 (0.03)
Epinephelus hexagonatus Hao 21.337 (0.907) 0.636 (0.068) -0.067 (0.043)
Epinephelus merra Moorea 14.078 (0.864) 1.192 (0.098) -0.016 (0.014)
Epinephelus merra Gambiers 18.52 (0.887) 0.655 (0.08) -0.093 (0.052)
Epinephelus merra Hao 22.481 (0.617) 0.842 (0.084) -0.043 (0.034)
Epinephelus polyphekadion Hao 38.812 (2.245) 0.248 (0.024) -0.332 (0.083)
Gnathodentex aureolineatus Gambiers 20.151 (0.881) 0.547 (0.076) -0.155 (0.079)
Halichoeres trimaculatus Gambiers 17.844 (1.584) 0.655 (0.157) -0.036 (0.058)
Lutjanus fulvus Moorea 21.82 (1.598) 1.175 (0.244) -0.019 (0.025)
Lutjanus fulvus Tuamotu 24.765 (1.386) 0.476 (0.076) -0.215 (0.11)
Mulloidichthys flavolineatus Moorea 25.762 (0.535) 1.601 (0.108) -0.008 (0.009)
Myripristis berndti Gambiers 23.324 (1.519) 0.338 (0.058) -0.391 (0.157)
Myripristis berndti Marquesas 21.795 (0.496) 0.378 (0.036) -0.298 (0.101)
Myripristis berndti Moorea 17.458 (0.9) 0.553 (0.064) -0.11 (0.044)
Naso lituratus Gambiers 31.83 (1.758) 1.718 (0.227) -0.005 (0.011)
Odonus niger Moorea 22.957 (1.147) 0.588 (0.071) -0.08 (0.041)
Ostorhinchus angustatus Moorea 7.216 (0.356) 0.686 (0.089) -0.11 (0.036)
Ostorhinchus apogonoides Marquesas 8.226 (0.217) 0.839 (0.053) -0.073 (0.014)
Plectropomus laevis Gambiers 77.388 (3.133) 0.253 (0.023) -0.408 (0.097)
Plectropomus laevis Hao 52.309 (2.535) 0.682 (0.054) -0.012 (0.019)
Pristiapogon taeniopterus Marquesas 8.025 (0.275) 0.778 (0.069) -0.093 (0.027)
Sargocentron microstoma Moorea 16.263 (0.542) 0.656 (0.052) -0.056 (0.027)
Siganus argenteus Gambiers 26.724 (2.244) 0.53 (0.098) -0.091 (0.067)
Stegastes nigricans Moorea 12.006 (0.886) 0.409 (0.099) -0.199 (0.11)

2.3 Growth curves (sp.)

# 1. Growth curves ----

read.csv("./../data/03_back-calculated_vbgf_fitted_sp.csv") %>% 
  ggplot(data = .) +
    geom_ribbon(aes(x = age, ymin = ypred_lq, ymax = ypred_uq), alpha = 0.5, fill = col_fill_graph, alpha = 0.6) +
    geom_line(aes(x = age, y = ypred_m), color = col_color_graph) +
    geom_point(data = data_complete %>%
                 filter(Species %in% unique(read.csv("./../data/03_back-calculated_vbgf_fitted_sp.csv")$Species)) %>% 
                 mutate(Li_sp_m = Li_sp_m/10), 
               aes(x = Agei, y = Li_sp_m), color = col_color_graph, fill = col_fill_graph, size = 1, shape = 21) +
    facet_wrap(~Species, scales = "free", ncol = 5) +
    lims(y = c(0, NA)) +
    theme(strip.text.x = element_text(face = "italic")) +
    labs(x = "Age (years)", y = "Length (TL, cm)") +
    guides(colour = guide_legend(override.aes = list(alpha = 1, size = 0.5)))
Figure 3. Growth curves on back-calculated data by species.

Figure 3. Growth curves on back-calculated data by species.

2.4 Growth curves (sp. loc.)

# 1. Growth curves ----

read.csv("./../data/03_back-calculated_vbgf_fitted_sploc.csv") %>% 
  ggplot(data = .) +
    geom_ribbon(aes(x = age, ymin = ypred_lq, ymax = ypred_uq, fill = Location), alpha = 0.5, show.legend = FALSE) +
    geom_line(aes(x = age, y = ypred_m, color = Location), show.legend = FALSE) +
    geom_point(data = data_complete %>%
                 filter(Species %in% unique(read.csv("./../data/03_back-calculated_vbgf_fitted_sploc.csv")$Species)) %>% 
                 mutate(Li_sploc_m = Li_sploc_m/10), 
               aes(x = Agei, y = Li_sploc_m, color = Location), size = 1, show.legend = FALSE) +
    facet_grid(Species~Location, scales = "free") +
    lims(y = c(0, NA)) +
    labs(x = "Age (years)", y = "Length (TL, cm)") +
    theme(strip.text.x = element_text(size = 8),
          strip.background = element_rect(colour = "black", fill = col_facet, size = 1),
          plot.title = element_text(colour = col_color_graph),
          plot.subtitle = element_text(colour = "black"),
          strip.text.y = element_text(angle = 360, face = "italic"))
Figure 4. Growth curves on back-calculated data by species and location.

Figure 4. Growth curves on back-calculated data by species and location.

3. Raw vs Back-calculation

3.1 Growth curves (sp.)

# 1. List of species with enough individuals for one location ----

species_list <- c("Abudefduf sexfasciatus", "Acanthurus achilles", "Acanthurus pyroferus", 
                  "Acanthurus triostegus", "Balistapus undulatus", "Centropyge flavissima",
                  "Chaetodon ornatissimus", "Chromis iomelas", "Ctenochaetus marginatus",
                  "Dascyllus aruanus", "Epinephelus hexagonatus", "Myripristis berndti",
                  "Odonus niger", "Ostorhinchus angustatus", "Ostorhinchus apogonoides",
                  "Pristiapogon taeniopterus", "Sargocentron microstoma", "Stegastes nigricans")

# 1. Transform the raw data ----

# 1.1 Points --

data_raw_points <- read.csv("./../data/02_back-calculated-size-at-age_morat-et-al.csv") %>%
  filter(Species %in% species_list) %>% 
  select(Family, Genus, Species, ID, Agecpt, Lcpt, Location, Observer) %>% 
  unique() %>%
  dplyr::group_by(Species, Location) %>%
  dplyr::mutate(n = length(unique(ID))) %>%
  filter(n >= 10) %>% # filter with at least 10 replicates
  ungroup() %>% 
  dplyr::mutate(Lcpt = Lcpt/10) %>%  # Convert to cm
  mutate(Type = "Population-level data")

# 1.2 Curves -- 

data_raw_curve <- read.csv("./../data/03_raw_vbgf_fitted_sp.csv") %>% 
  filter(Species %in% species_list) %>% 
  mutate(Type = "Population-level data")

# 2. Transform the back-calculated data ----

# 2.1 Points --

data_backcalc_points <- read.csv("./../data/02_back-calculated-size-at-age_morat-et-al.csv") %>% 
  filter(Species %in% species_list) %>% 
  dplyr::mutate(Li_sp_m = Li_sp_m/10) %>%  # Convert to cm
  mutate(Type = "Back-calculated data")

# 2.2 Curves --

data_backcalc_curve <- read.csv("./../data/03_back-calculated_vbgf_fitted_sp.csv") %>% 
  filter(Species %in% species_list) %>% 
  mutate(Type = "Back-calculated data")

# 2. Make the plot ----

ggplot() +
  # Raw data
  geom_ribbon(data = data_raw_curve, aes(x = Agei, ymin = Q2.5, ymax = Q97.5, fill = Type), alpha = 0.5) +
  geom_line(data = data_raw_curve, aes(x = Agei, y = Estimate, color = Type), show.legend = FALSE) +
  geom_point(data = data_raw_points, aes(x = Agecpt, y = Lcpt, color = Type), size = 1) +
  # Back-calculated data
  geom_ribbon(data = data_backcalc_curve, aes(x = age, ymin = ypred_lq, ymax = ypred_uq, fill = Type), alpha = 0.5) +
  geom_line(data = data_backcalc_curve, aes(x = age, y = ypred_m, color = Type), show.legend = FALSE) +
  geom_point(data = data_backcalc_points, aes(x = Agei, y = Li_sp_m, color = Type), size = 1) +
  # Appearance
  facet_wrap(~Species, scales = "free", ncol = 3) +
  lims(y = c(0, NA)) +
  theme(strip.text.x = element_text(face = "italic"),
        legend.position = "top",
        legend.title = element_blank()) +
  labs(x = "Age (years)", y = "Length (TL, cm)") +
  lims(y = c(0, NA), x = c(0, NA)) +
  scale_color_manual(values = c(col_color_graph, "#d91e18")) +
  scale_fill_manual(values = c(col_fill_graph, "#e74c3c")) +
  guides(colour = guide_legend(override.aes = list(alpha = 1, size = 0.5)))
Figure 5. Comparison of growth curves by species between back-calculated data and population-level data.

Figure 5. Comparison of growth curves by species between back-calculated data and population-level data. 

# 3. Remove useless objects ----

rm(data_raw_points, data_raw_curve, data_backcalc_points, data_backcalc_curve)

3.2 Growth curves (sp. loc.)

# 1. Transform the raw data ----

species_list2 <- species_list[!species_list %in% c("Balistapus undulatus", "Centropyge flavissima", "Chaetodon ornatissimus")]

# 1.1 Points --

data_raw_points <- read.csv("./../data/02_back-calculated-size-at-age_morat-et-al.csv") %>%
  filter(Species %in% species_list2) %>% 
  select(Family, Genus, Species, ID, Agecpt, Lcpt, Location, Observer) %>% 
  unique() %>%
  dplyr::group_by(Species, Location) %>%
  dplyr::mutate(n = length(unique(ID))) %>%
  filter(n >= 10) %>% # filter with at least 10 replicates
  ungroup() %>% 
  dplyr::mutate(Lcpt = Lcpt/10) %>%  # Convert to cm
  mutate(Type = "Raw")

# 1.2 Curves -- 

data_raw_curve <- read.csv("./../data/03_raw_vbgf_fitted_sploc.csv") %>% 
  filter(Species %in% species_list2) %>% 
  mutate(Type = "Raw")

# 2. Transform the back-calculated data ----

# 2.1 Points --

data_backcalc_points <- read.csv("./../data/02_back-calculated-size-at-age_morat-et-al.csv") %>% 
  filter(Species %in% species_list2,
         !(Species == "Abudefduf sexfasciatus" & Location == "Moorea"),
         !(Species == "Acanthurus triostegus" & Location == "Marquesas"),
         !(Species == "Myripristis berndti" & Location == "Gambiers"),
         !(Species == "Myripristis berndti" & Location == "Marquesas")) %>% 
  dplyr::mutate(Li_sploc_m = Li_sploc_m/10) %>%  # Convert to cm
  mutate(Type = "Back-calculated")

# 2.2 Curves --

data_backcalc_curve <- read.csv("./../data/03_back-calculated_vbgf_fitted_sploc.csv") %>% 
  filter(Species %in% species_list2,
         !(Species == "Abudefduf sexfasciatus" & Location == "Moorea"),
         !(Species == "Acanthurus triostegus" & Location == "Marquesas"),
         !(Species == "Myripristis berndti" & Location == "Gambiers"),
         !(Species == "Myripristis berndti" & Location == "Marquesas")) %>% 
  mutate(Type = "Back-calculated")

# 2. Make the plot ----

ggplot() +
  # Raw data
  geom_ribbon(data = data_raw_curve, aes(x = Agei, ymin = Q2.5, ymax = Q97.5, fill = Type), alpha = 0.5) +
  geom_line(data = data_raw_curve, aes(x = Agei, y = Estimate, color = Type), show.legend = FALSE) +
  #geom_point(data = data_raw_points, aes(x = Agecpt, y = Lcpt, color = Type), size = 1) +
  # Back-calculated data
  geom_ribbon(data = data_backcalc_curve, aes(x = age, ymin = ypred_lq, ymax = ypred_uq, fill = Type), alpha = 0.5) +
  geom_line(data = data_backcalc_curve, aes(x = age, y = ypred_m, color = Type), show.legend = FALSE) +
  #geom_point(data = data_backcalc_points, aes(x = Agei, y = Li_sploc_m, color = Type), size = 1) +
  # Appearance
  facet_wrap(~Species, scales = "free", ncol = 3) +
  lims(y = c(0, NA)) +
  theme(strip.text.x = element_text(face = "italic"),
        legend.position = "top",
        legend.title = element_blank()) +  labs(x = "Age (years)", y = "Length (TL, cm)") +
  lims(y = c(0, NA), x = c(0, NA)) +
  scale_color_manual(values = c(col_color_graph, "#d91e18")) +
  scale_fill_manual(values = c(col_fill_graph, "#e74c3c")) +
  guides(colour = guide_legend(override.aes = list(alpha = 1, size = 0.5)))
Figure 6. Comparison of growth curves by species and location (one location only by growth curve) between back-calculated data and population-level data.

Figure 6. Comparison of growth curves by species and location (one location only by growth curve) between back-calculated data and population-level data. 

# 3. Save the plot ----

ggsave("./../figs/figure-2_comparison-growth-curves-raw-vas-back-calculated.png", height = 15, width = 10) # PNG
ggsave("./../figs/figure-2_comparison-growth-curves-raw-vas-back-calculated.eps", height = 15, width = 10, device = cairo_pdf) # PDF

# 4. Remove useless objects ----

rm(data_raw_points, data_raw_curve, data_backcalc_points, data_backcalc_curve)

4. Comparison with literature

# 1. Import file synthesizing the VBGF parameters from literature ----

vbgf_literature <- read_excel("./../data/00_von-bertalanffy-literature.xlsx", sheet = 1) %>% 
  dplyr::mutate(Linf = ifelse(Size_unit == "mm", Linf/10, Linf),
                Size_max = ifelse(Size_unit == "mm", Size_max/10, Size_max)) %>% # Convert all length values to cm
  select(-Family, -Genus) # Remove Family and Genus to add those level through fishbase

# 2. Check validity of species names ----

vbgf_literature <- vbgf_literature %>% 
  rowwise() %>% 
  dplyr::mutate(Species2 = ifelse(is_empty(validate_names(Species)) == TRUE, 
                           NA, 
                           validate_names(Species))) %>% # Correct the names through fishbase
  ungroup() %>% 
  dplyr::mutate(Species = str_replace_all(Species, c("Acanthurus chirugus" = "Acanthurus chirurgus",
                                                     "Scarus psitticus" = "Scarus psittacus"))) %>% # Correct unfound names
  filter(Species %in% unique(data_complete$Species)) %>% 
  mutate(Reference = paste0(Reference, " (", Reference_type, ")")) %>%
  select(-Size_unit, -Sex, -Reference_type, -Species2) %>% 
  arrange(Species) %>% 
  filter(Species %in% unique(read.csv("./../data/03_back-calculated_vbgf_predictions_sp.csv")$Species))

4.1 Growth parameters

4.1.1 Table

# 1. Get age max and length max by species ----

table_agemax <- data_complete %>% 
  drop_na(Li_sp_m) %>% 
  select(Species, Agecpt, Lcpt) %>% 
  group_by(Species) %>% 
  summarise(Age_max = max(Agecpt, na.rm = TRUE),
            Size_max = round(max(Lcpt, na.rm = TRUE)/10, 3))

# 2. Join mean estimated parameters, with size and length max, and with literature ----

read.csv("./../data/03_back-calculated_vbgf_predictions_sp.csv") %>% 
  mutate(Estimate = round(Estimate, 3)) %>% 
  pivot_wider(c("Species"), names_from = Parameter, values_from = Estimate) %>% 
  select(Species, linf, k, t0) %>% 
  rename(Linf = linf, K = k) %>% 
  mutate(Reference = "This study") %>% 
  left_join(., table_agemax) %>% 
  full_join(., vbgf_literature) %>%
  select(Species, Linf, K, t0, Age_max, Size_max, Reference) %>% 
  arrange(Species) %>% 
  write.csv2(., "./../figs/online-only-table-2_comparison-vbgf-parameters-literature.csv", row.names = FALSE)

4.1.2 Plot

# 1. Comparison of VBGF parameters between our study and literature ----

vbgf_literature %>% 
  select(Species, Linf, K, t0, Reference) %>% 
  pivot_longer(c("Linf", "K", "t0"), names_to = "Parameter", values_to = "Estimate") %>% 
  mutate(Estimate = as.numeric(Estimate),
         Type = "Literature") %>% 
  bind_rows(., read.csv("./../data/03_back-calculated_vbgf_predictions_sp.csv") %>% 
                select(-Est.Error, -Q2.5, -Q97.5) %>% 
                mutate(Parameter = str_replace_all(Parameter, c("linf" = "Linf",
                                                                "k" = "K")),
                       Type = "This study",
                       Reference = "This study")) %>% 
  mutate(Species = as.factor(Species),
         Parameter = factor(Parameter, levels = c("Linf", "K", "t0"))) %>% 
  ggplot(data = ., aes(x = reorder(Species, desc(Species)), y = Estimate, color = Type)) +
    geom_point() +
    facet_wrap(~Parameter, scales = "free") +
    coord_flip() +
    labs(x = NULL, y = "Value") +
    theme(axis.text.y = element_text(face = "italic"))
Figure 8. Comparison of VBGF growth parameters estimated through Bayesian framework on species back-calculated data, with parameters from literature.

Figure 8. Comparison of VBGF growth parameters estimated through Bayesian framework on species back-calculated data, with parameters from literature.

4.2 Growth curves

# 1. Get predictions from literature data ----

vbgf_literature_pred <- vbgf_literature %>% 
  # Complete missing Age_max
  filter(is.na(Age_max)) %>% 
  left_join(., vbgf_literature %>% 
              group_by(Species) %>% 
              summarise(Age_max2 = max(Age_max, na.rm = TRUE)) %>% 
              mutate(Age_max2 = na_if(Age_max2, -Inf))) %>% 
  mutate(Age_max = ifelse(is.na(Age_max), Age_max2, Age_max)) %>% 
  select(-Age_max2) %>% 
  bind_rows(., vbgf_literature %>% 
              filter(!(is.na(Age_max)))) %>% 
  arrange(Species) %>% 
  # Fix missing t0 to 0
  dplyr::mutate(t0 = ifelse(is.na(t0), 0, t0)) %>% 
  filter(!(is.na(K)), !(is.na(Linf)), !(is.na(t0))) %>% 
  mutate_at(c("Linf", "K", "t0"), as.numeric) %>% 
  filter(Species %in% unique(read.csv("./../data/03_back-calculated_vbgf_fitted_sp.csv")$Species)) %>% 
  # Make the predictions
  dplyr::mutate(Line = row_number()) %>% 
  group_by(Line) %>% 
  do(pred_vbgf(data = .))

# 2. Make the plot ----

plot_final <- read.csv("./../data/03_back-calculated_vbgf_fitted_sp.csv") %>% 
  filter(Species %in% unique(vbgf_literature_pred$Species)) %>% 
  ggplot(data = .) +
  geom_line(data = vbgf_literature_pred, aes(x = Agei, y = Li, group = as.factor(Line), 
                                             text = Reference), color = "#6c7a89", size = 0.7) +
  geom_ribbon(aes(x = age, ymin = ypred_lq, ymax = ypred_uq), alpha = 0.5, fill = col_fill_graph) +
  geom_line(aes(x = age, y = ypred_m), size = 0.7, color = col_color_graph) +

  facet_wrap(~Species, ncol = 4, scales = "free") +
  labs(x = "Age (years)", y = "Length (cm)") +
  theme(strip.text.x = element_text(face = "italic")) +
  lims(y = c(0, NA), x = c(0, NA)) +
  guides(colour = guide_legend(override.aes = list(alpha = 1, size = 0.5)))

plot_final
Figure 9. Comparison of growth curves based on VBGF growth parameters estimated through Bayesian framework on species back-calculated data with growth curves from literature.

Figure 9. Comparison of growth curves based on VBGF growth parameters estimated through Bayesian framework on species back-calculated data with growth curves from literature. 

# 4. Convert to plotly ----

plot_final <- plot_final +
  theme(legend.position = "none")

ggplotly(plot_final, tooltip = c("Reference", "Location", "Size_type"))

Figure 10. Comparison of growth curves based on VBGF growth parameters estimated through Bayesian framework on species back-calculated data with growth curves from literature (interactive plot).

5. Annex

# 1. Merge the tables on growth parameters by species and by species and location in one table ----

table_final <- read.csv("./../data/03_back-calculated_vbgf_predictions_sp.csv") %>% 
  mutate(Location = "All") %>% 
  bind_rows(., read.csv("./../data/03_back-calculated_vbgf_predictions_sploc.csv")) %>% 
  select(-Est.Error) %>% 
  pivot_wider(names_from = Parameter, names_glue = "{Parameter}_{.value}", values_from = c(Estimate, Q2.5, Q97.5)) %>% 
  arrange(Species, Location) %>% 
  mutate_if(is.numeric, ~round(.x, digits = 3)) %>% 
  mutate(linf_range = paste0(linf_Q2.5, " ; ", linf_Q97.5),
         k_range = paste0(k_Q2.5, " ; ", k_Q97.5),
         t0_range = paste0(t0_Q2.5, " ; ", t0_Q97.5)) %>% 
  select(Species, Location, linf_Estimate, linf_range, k_Estimate, k_range, t0_Estimate, t0_range) %>% 
  # Remove rows of only one sploc (-> sp sufficient) 
  group_by(Species) %>% 
  mutate(n = n()) %>% 
  filter(n > 2 | Location == "All") %>% 
  select(-n)

# 2. Create the table of number of individuals used to estimates parameters ----

# 2.1 Table for sp --

table_n_sp <- data_complete %>% 
  drop_na(Li_sp_m) %>%
  group_by(ID) %>%
  filter(Agecpt > 2) %>% # filter with maxage > 2
  ungroup() %>%
  dplyr::group_by(Species) %>%
  dplyr::summarise(n = length(unique(ID))) %>%
  filter(n >= 5) %>% # filter with at least 5 individuals
  ungroup() %>% 
  mutate(Location = "All")

# 2.1 Table for sploc --

table_n_sploc <- data_complete %>% 
  drop_na(Li_sploc_m) %>%
  group_by(ID) %>%
  filter(Agecpt > 2) %>% # filter with maxage > 2
  filter(!(ID %in% c("SA_MI_MA_03_17_128", "CE_FL_MA_03_17_214", "CE_FL_MA_03_17_216"))) %>% # Filter ID leading to non-convergence
  ungroup() %>%
  dplyr::group_by(Species, Location) %>%
  dplyr::summarise(n = length(unique(ID))) %>%
  filter(n >= 5) %>% # filter with at least 5 individuals
  ungroup()

# 2.3 Merge the n tables --

table_n <- full_join(table_n_sp, table_n_sploc) %>% 
  arrange(Species, Location)

# 3. Merge the n table with the 

table_final %>% 
  left_join(., table_n) %>% 
  select(Species, Location, n, linf_Estimate, linf_range, k_Estimate, k_range, t0_Estimate, t0_range) %>% 
  write.csv2(., 
            "./../figs/online-only-table-2_von-bertalanffy-parameters-estimates.csv", 
            row.names = FALSE)

Reproducibility

# 1. Reproducibility ----

sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
## [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
## [5] LC_TIME=French_France.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rfishbase_3.0.4     plotly_4.9.2.1      readxl_1.3.1       
##  [4] formattable_0.2.0.1 kableExtra_1.1.0    forcats_0.5.0      
##  [7] stringr_1.4.0       dplyr_1.0.0         purrr_0.3.4        
## [10] readr_1.3.1         tidyr_1.1.0         tibble_3.0.1       
## [13] ggplot2_3.3.2       tidyverse_1.3.0     extrafont_0.17     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.4.6      lubridate_1.7.9   lattice_0.20-41   assertthat_0.2.1 
##  [5] digest_0.6.25     R6_2.4.1          cellranger_1.1.0  backports_1.1.7  
##  [9] reprex_0.3.0      evaluate_0.14     highr_0.8         httr_1.4.1       
## [13] pillar_1.4.4      rlang_0.4.6       curl_4.3          lazyeval_0.2.2   
## [17] data.table_1.12.8 rstudioapi_0.11   extrafontdb_1.0   blob_1.2.1       
## [21] rmarkdown_2.3     labeling_0.3      webshot_0.5.2     htmlwidgets_1.5.1
## [25] munsell_0.5.0     broom_0.5.6       compiler_4.0.1    modelr_0.1.8     
## [29] xfun_0.14         pkgconfig_2.0.3   htmltools_0.5.0   tidyselect_1.1.0 
## [33] fansi_0.4.1       viridisLite_0.3.0 crayon_1.3.4      dbplyr_1.4.4     
## [37] withr_2.2.0       grid_4.0.1        nlme_3.1-148      jsonlite_1.6.1   
## [41] Rttf2pt1_1.3.8    gtable_0.3.0      lifecycle_0.2.0   DBI_1.1.0        
## [45] magrittr_1.5      scales_1.1.1      cli_2.0.2         stringi_1.4.6    
## [49] farver_2.0.3      fs_1.4.1          xml2_1.3.2        ellipsis_0.3.1   
## [53] generics_0.0.2    vctrs_0.3.1       gh_1.1.0          tools_4.0.1      
## [57] glue_1.4.1        crosstalk_1.1.0.1 hms_0.5.3         yaml_2.2.1       
## [61] colorspace_1.4-1  rvest_0.3.6       memoise_1.1.0     knitr_1.28       
## [65] haven_2.3.1

Jeremy WICQUART | | 2020-08-21 14:06:31